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Abstract. We test through stellar N-body simulations 
some scenarios to explain the dynamics of the peculiar 
nucleus of the Andromeda galaxy (M 31): although HST 
observations reveal a double nucleus morphology, the ro¬ 
tation field is almost symmetric around the bulge gravity 
centre and the velocity dispersion is off-centred. We show 
that any to = 1 perturbation has a very short life-time (a 
few 10® yr). Assuming that the bright peak (PI) is a cold 
stellar cluster infalling into the nucleus, and that the large 
central velocity gradient is due to a central dark mass (in 
the range 7 10^-10® Mq), we obtain a reasonably good 
fit to the observations. However, if this cluster lies in the 
central 20 pc, we estimate the life-time of the cluster to be 
less than 0.5 Myr. The dynamical friction is more efficient 
than estimated by analytic formulae, and is essentially due 
to the deformation of the stellar cluster through the huge 
tidal forces provided by the black hole. We show that the 
cluster cannot be on a circular orbit around the centre if 
the nucleus hosts a massive black hole of a few lO’’ Mq, 
and finally provide some estimates of the kinematics as 
observed with HST. 

Key words: galaxies: M 31 - galaxies: nuclei - galaxies: 
kinematic and dynamics 


1. Introduction 

Each new observational result concerning the nucleus of 
M 31 seems to deepen the mystery of its structure. This 
is mainly due to the improvement in the achieved spa¬ 
tial resolution of the photometric and spectroscopic data: 
the better resolved the object is, the more complex it ap¬ 
pears. For any proposed physical mechanism which can 
quantitatively reproduce the observables, it is important 
to provide some specific predictions. 

Send offprint requests to: E. Emsellem (email: 
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In this respect, the nucleus of M 31 is certainly an ex¬ 
cellent laboratory to test our knowledge of galactic nuclei, 
as we can now optically resolve scales as small as 0.3 pc 
(e.g. HST or optically adaptive ground-based systems). 
Indeed, the numerous studies achieved on this object re¬ 
vealed a puzzling complexity. 

M 31’s nucleus is a compact stellar system with an 
average ellipticity of (e) ^ 0.4. Asymmetries in its major- 
axis surface brightness profile (in the visible) were already 
observed in 1974 by Light et ah, and subsequently under¬ 
stood using the pre COSTAR HST data as a double com¬ 
ponent nucleus (Lauer et al. 1993). In the V band, the 
fainter peak (P2) is almost coincident with the centre of 
the bulge isophotes, as the brightest one (PI) is located 
at about 0/5 (~ 1.8 pc for a distance to M 31 of 0.77 
Mpc) from P2 roughly along the major-axis. At 175 nm 
however, P2 is the brightest point, with PI having an UV 
upturn similar to its surroundings (King et al. 1995). 

In 1960, Lallemand et al. obtained spectrograms of the 
central region of M 31 and detected the rapid rotation of 
its nucleus, concluding that it is a dynamically indepen¬ 
dent structure. This was confirmed by Kormendy (1988) 
and Dressier & Richstone (1988) by long-slit spectroscopy 
with CCD detectors. Both studies suggest the presence of 
a central dark mass of the order of 10^ Mq in the centre 
of the nucleus which would account for the high value of 
the central velocity and velocity dispersion gradients. 

Sub-arcsecond velocity and dispersion maps were ob¬ 
tained with the TIGER spectrograph by Bacon et al. 
(1994, hereafter BEMN94) which uncovered more asym¬ 
metries (already traced in previous published data): the 
velocity field is nearly symmetric about a point Vq located 
very close to P2 (0/05 according to Lauer et al. 1993), but 
the maximum S of the stellar velocity dispersion roughly 
corresponds to the symmetric point of PI with respect to 
P2 (see Fig. 26 of BEMN94). 

Different interpretations have been proposed to explain 
the observed photometric and dynamical asymmetries in 
the nucleus of M 31 (see Sect. ^ and references therein), 
but they each encounter severe problems. Some questions 
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remain open, such as the actual efficiency of dynamical 
friction for an external body falling into the nucleus, the 
life-time of m = 1 perturbations, or the amount of dynam¬ 
ical perturbation on the nucleus itself; we try here to an¬ 
swer these questions through N body simulations. In this 
paper, we concentrate on one of the proposed mechanism, 
namely that PI corresponds to a stellar cluster falling into 
the potential well of the nucleus, hosting a massive black 
hole. In a companion paper, we will examine whether a 
central black hole is unavoidable, or the presence of a nu¬ 
clear bar could account for the observational data. 

We summarize and comment the alternatives already 
discussed by different authors in Sect. 2. In Sect. 3 we de¬ 
scribe the N body code we used as well as the initial con¬ 
ditions of the different experiments. The results are dis¬ 
cussed in Sect. 4. In Sect. 5, we present and discuss our 
attempts to fit the observables using different assump¬ 
tions. The general discussion and conclusions are given 
in Sect. 6. 

2. Different mechanisms 

The M31 nucleus is an object of peculiar attention, be¬ 
cause of its large central velocity gradient and velocity 
dispersion: it is one of the best candidate among nearby 
galaxies to host a massive black hole (BH of the order 
of 10^-10® Mq, e.g. Kormendy & Richstone 1995). Possi¬ 
ble scenarii include one massive BH, two BHs, or no BH 
at all. We will consider here only the first class of mod¬ 
els, with several possibilities to account for the observed 
asymmetries. Previous studies already discarded a grav¬ 
itational lens or a dust component to be responsible for 
these asymmetries (e.g. BEMN94), so we will discuss the 
two remaining ones: a dense stellar cluster falling into the 
nucleus (already suggested by Dressier & Richstone 1988, 
see also Lauer et al. 1993, BEMN94) or a thick eccentric 
disc of stars on Keplerian orbits around the black hole, 
maintained through self-gravity and dynamical friction on 
the bulge (Tremaine 1995). 

2.1. Is the disc eccentric ? 

Tremaine (1995) recently proposed a model based on a 
set of aligned Keplerian orbits, forming an m = 1 pattern. 
The main argument for an orbit crowding effect is the lack 
of any significant colour gradient between PI and the rest 
of the nucleus, except for P2 which has a higher UV up¬ 
turn which could correspond to the radio source detected 
by Crane et al. (1992) and linked with the presumed cen¬ 
tral dark mass. But as argued by King et al. (1995), this 
difference in the UV flux is equivalent to a single PAGB 
star, although it seems more extended (P. Crane, priv. 
comm.). Tremaine (1995) also argued that the observed 
isophote twist of PI as well as the elongations of PI and 
P2 along the P1-P2 axis are natural consequences of the 
eccentric disc hypothesis. The only problem in Tremaine’s 


model lies in the strong predicted asymmetry of the ob¬ 
served velocity field, in contradiction with the recently 
acquired data of Kormendy (SIS/CFHT, ct* ^ 0'.'27). 

In his proposition, Tremaine (1995) suggests that the 
alignment of apsides of the orbits forming the eccentric 
disc is maintained by the self-gravity of the nuclear disc. 
In fact, an m = 1 mode could amplify in the disc, if it 
has a positive pattern speed Up, and if dynamical friction 
against the non-rotating bulge removes effectively the an¬ 
gular momentum, provoquing more eccentricity in the disc 
orbits. The efficiency of this mechanism is still unknown. 

2.2. An m = 1 density wave 

Another possibility is that the m = 1 perturbation is tran¬ 
sient, but slowly damped. Weinberg (1994) shows that a 
stellar system can sustain weakly damped m = 1 mode 
for hundreds of crossing times. A fly-by encounter could 
excite such a mode, and explain off-centring in some spiral 
galaxies. We have tried to estimate the relevant time-scale 
for the fading of an m = 1 excitation in the particular con¬ 
ditions of the M31 nucleus, in the presence of a massive 
black hole. The excitation was provided by a slight initial 
displacement of the black hole, either in space or velocity. 
An m = 1 spiral wave was generated in the nuclear disc 
(see Fig. 1^, which disappeared through phase mixing in 
about 0.5 Myr, while the black hole was braked down to 
the centre by dynamical friction. The perturbation kept 
a positive pattern speed Q= 50km/s/pc during about 3-4 
rotations (the corresponding period is 0.125 Myr). After 
that only a material arm was precessing at the negative 
velocity -3.6km/s/kpc until 1 Myr. Finally, an unsta- 



-2.0 0-0 2.0 -2.0 0.0 2.0 
X (arcsec) x (arcsec) 


Fig. 1. Projected isophotes of the nuclear disc with the black 
hole launched with an initial offset. Left:face-on, Right:edge-on; 
Top: 0.08 Myr; Bottom: 0.32 Myr. 

ble m = 1 mode may also exist in the presence of gas in 
the centre of M31 (Shu et al, 1990, Junqueira & Combes 
1996). 







E. Emsellem & F. Combes: N body simulations of the nucleus of M 31 


3 


2.3. A falling stellar cluster 

Even at HST resolution, the nucleus of M 31 appears 
as the superposition of two rather smooth components 
corresponding to PI and P2. The detailed structure of 
each individual components depends on the decomposi¬ 
tion method (BEMN94, King et al. 1995) but the corre¬ 
sponding global parameters seem to vary very little (see 
Table 1^). Let us emphasize, as already did King et al. 
(1995), that although the brightest point of the nuclear 
region of M 31 corresponds to PI, it never contributes 
more than 0.55 of the total surface brightness (assuming 
P2 is symmetric) and the integrated magnitude of PI is 
2.6 magnitude fainter than of P2. Occam’s razor suggests 


Table 1. Global parameters for PI and P2: absolute magni¬ 
tude in the V band, core radius, mean ellipticity, relative po¬ 
sition angle (PA in degrees), and displacement (in arcsecond) 
along the nucleus major {Sx') and minor-axis {Sy'). 


# 

Mv 

Vc 

(^> 

(PA) 

5x' 

5y' 


(mag) 

(") 


(“) 

(") 

(") 

PI 

-9.6 

0.25 

0.31 

104 

-0.5 

0.14 

P2 

-12.2 

0.88 

0.39 

54 

0 

0 


to favor the hypothesis of P2 being a nearly symmetric 
system, and PI a superimposed elongated cluster. Since 
the velocity profile does not seem to be highly perturbed 
by the presence of PI (even at the resolution of u* ~ 0'.'27 
attained by SIS data), it suggests that the line-of-sight ve¬ 
locity of PI is close to the line-of-sight velocity of the sur¬ 
rounding nuclear stars (including the seeing effect). In this 
scenario, the main effect is due to the low dispersion of the 
superimposed cluster which produces the offset observed 
in the velocity dispersion field (BEMN94). We examine 
this point in details in Sects || and ||. 

3. N-body models 

3.1. The code 

We have simulated the infalling of a compact stellar 
system in the nucleus of M 31 using an N-body code 
whose main characteristics were described in Combes et al. 
(1990). A 3D cartesian grid was used, with free boundary 
conditions, and we used the method described by James 
(1977) to avoid a 2^ multiplication of the number of cells 
in the grid. The total number of particles was 152384, and 
the grid size 128 x 128 x 80. The time step of integration 
was 200 yr. 

3.2. Initial conditions 

As already suggested by Lallemand (1960), the nucleus 
of M 31 appears as a distinct photometric and dynami¬ 


cal component. Indeed, observed properties of the galaxy 
exhibit strong changes as we reach a radius of about 3 
arcseconds along the major-axis: 

— the surface brightness rises by more than 1.5 magni¬ 
tude in the central 2 arcseconds, while the outer major- 
axis profile of the bulge only shows an increase of 0.1 
magnitude in the same interval. 

— the stellar velocity is nearly zero at 3 arcseconds but 
peaks around 155 km.s~^at 078 (cr* ~ 0727, Kormendy 
& Richstone 1995, hereafter KR95). 

— the stellar velocity dispersion stays nearly constant in 
the inner bulge with cr ^ 150 km.s”^, and starts to 
rise in the inner 2 arcseconds to reach a value of ^ 250 
km.s~^(KR95) in the centre. 

Moreover, Bacon et al. (1994) have shown that the 
assumption of the nucleus being a dynamically isolated 
entity leads to differences of less than 1% in the deriva¬ 
tion of a local tensor virial theorem. This encouraged us 
to treat the bulge and the nucleus as two separate com- 
ponenta.. We therefore adapted the grid to the size of the 
nucleusEJ, by taking a cell size of 0.165 pc which led to a 
volume of about 21 x 21 x 13 pc^ centred on P2. 

3.2.1. The bulge 

Since the nucleus appears as a dynamically isolated com¬ 
ponent, we include the contribution of the bulge as a fixed 
potentialEl. 

We have assumed the bulge to be spherically symmet¬ 
ric in our simulations. This assumption leads to an error of 
less than 10% in the potential of the bulge inside the vol¬ 
ume of the grid. Moreover, in the presence of a supermas- 
sive black hole of a few 10^ Mq, this error is negligible in 
the central part of the nucleus. We thus approximated the 
mass density distribution of the bulge by a sum of Plum¬ 
mer spheres of different masses and scales: these functions 
yield simple analytic formula for the corresponding grav¬ 
itational potential and forces. 



where Mi and bt are the total mass and the scale of the 
component respectively. We have fitted the luminos¬ 
ity profile of the bulge using the axisymmetric model of 
BEMN94 (see Table 11) and the data of Kent (1983). This 
was achieved with a set of 4 Plummer spheres whose pa¬ 
rameters are given in Table H. Fig|| shows the excellent 

^ Throughout this paper, we use a distance of 0.77 Mpc for 
M 31 which leads to a scale of l" ~ 3.73 pc. 

^ Note that the definition of the “bulge” depends on the de¬ 
composition method. We rely here on the bulge/nucleus de¬ 
composition given in BEMN94. 
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agreement between this model and the published data 
on the inner bulge. The fact that we include the bulge 



Fig. 2. Major-axis surface brightness profile in the I band 
of the fixed bulge compared to the inner axisymmetric bulge 
model of Bacon et al. (1994, squares), and the values given by 
Kent outside 40" (circles, rescaled for the I band). 


Table 2. Parameters of the Plummer spheres for the bulge 
component: total mass {Mi), Plummer scale {bi) and the inte¬ 
grated mass inside a sphere of 10 pc {M}^). 


i 

Mi 

(Mo) 

bi 

(pc) 

(Mo) 

1 

3.7 10' 

16.5 

5.2 10® 

2 

3.0 10® 

55 

1.7 10® 

3 

3.3 lO'^ 

176 

6.0 10^ 

4 

3.7 10“ 

770 

8.1 10® 


contribution through a fixed potential allowed us to use 
all 152384 particles for the nucleus and the falling stel¬ 
lar cluster. Following BEMN94, we fixed the total visi¬ 
ble mass of both components to 1.9 10^ M©, leading to 
a mass of ~ 125 MQper particle. The mass of the bulge 
inside the grid can easily be derived with the formula: 
M{r) = Mir^/ (r^ -|- )^/^. Inside r = 10 pc the bulge 
contributes for a mass of ~ 7.6 10^ M!©: which is negligible 
compared to the added masses of PI and P2. 

3.2.2. The nucleus 

The nucleus is represented by 138530 particles (~ 91% of 
the total number), with an equivalent mass of ~ 1.73 10^ 
M©. We simulated it as a Toomre disc of size 3.8 pc. 
The vertical structure follows initially a sech^ distribution 
with a vertical scale of 1.7 pc at the origin, and slightly 


increasing with radius to fit the observed photometry. 
The radial variation of height with radius was chosen as 
(1 -I- (r/9pc)^)°'^^. The radial velocity dispersion is taken 
as an exponential law, to better approach the data: 

ar = (Joexp{-r/ry) 

with tTo=280km/s, and r„ = 5pc. The azimuthal disper¬ 
sion is taken, following the epicyclic approximation, as 
ae = Since the potential is nearly spherical to¬ 

wards the centre due to the presence of the massive black 
hole that dominates the mass, a first approximation of 
the z dispersion is to take it equal to the radial disper¬ 
sion az = (Jr- However, since this choice was not enough 
to keep, after relaxation, a nucleus as thick as observed, 
our final choice for the nearly edge-on nuclear disc was 
az = 1.3CTr. The actual rotational velocity was computed, 
after subtracting the asymmetric drift, according to Jeans 
equations (cf Binney & Tremaine 1987). 

3.2.3. The stellar cluster 

The stellar cluster is represented with the remaining 13854 
particles (9% of the particles corresponding to ^ 1.73 10® 
M©), as a truncated Plummer sphere with an equivalent 
core radius Cc = 0.13 pc (initial rc to anticipate any evolu¬ 
tion). The tidal radius Rt of the stellar cluster is adjusted 
to its initial environment, since it is launched at a radius 
i?o=6-9pc from the centre {Rt is between 1.6 and 2.2pc). 
The resulting concentrations c = log(i?t/rc) ~ 1 — 1.2, 
fall within the whole range of concentrations observed for 
globular clusters in our Galaxy (0.75 to 1.75). The veloci¬ 
ties follow an isotropic maxwellian distribution truncated 
at the escape velocity, and with a temperature parameter, 
providing the amount of kinetic energy corresponding to 
the virial theorem. 

The cluster is launched in the equatorial plane (at the 
radius i?o=6-9pc) with the velocity Vq (see Table ||). In 
most of the experiments, Vq corresponds to the circular 
velocity at Rq. This was motivated by two facts: first that 
dynamical friction does not depend much on the orbit his¬ 
tory of the satellite (e.g. Lin & Tremaine 1983) provided it 
is not yet disrupted. Since the central density of the clus¬ 
ter is an order of magnitude higher than 10® M© pc“®, 
the average density within 6pc including the black hole, 
we know that the tidal disruption has not occured at this 
radius (only the envelope of the cluster has been stripped 
away, which is taken into account through the initial Rt). 
And second, dynamical friction has the tendency to circu¬ 
larise orbits (Bontekoe & van Albada 1987). However, we 
have in Sect. ^ also performed experiments with a cluster 
on an hyperbolic orbit: this will give us insights for the 
particular case of M31’s nucleus (see Sect. ||). 

3.2.4. The black hole 

The observed steep central velocity gradient and veloc¬ 
ity dispersion peak suggested the presence of a central 
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dark mass of a few 10^ Mq (Dressier & Richstone 1988, 
Kormendy 1988). Bacon et al. (1994) confirmed this re¬ 
sult with self-gravitating two-integral dynamical models 
of the nucleus of M 31, and derived a central dark mass of 
~ 7.10^ 1^0 • We model it with a Plummer sphere with a 
core radius rc = 0.13 pc (i.e. limited by our resolution of 
0.15 pc) and a mass of a few times lO”^ MQ(see Table ||). 

3.3. Initial parameters 

We restrict our analysis to some of the parameters we 
thought to be crucial for the evolution of the system (cf 
Table I): 

— the black hole mass M^h, 

— the mass-to-light ratio of the cluster {M/L)c (i.e. its 
mass), 

— the initial radius of the cluster i?o- 

— the initial velocity of the cluster Vq. 


Table 3. Parameters of the different experiments: the black 
hole mass Mhh (Mq), the ratio between the M/L of the nu¬ 
cleus and the cluster ri{M/L), the core radius of the cluster Vc 
(pc), the radius Ro (pc) and velocity Vo (km.s“^) at which the 
cluster is launched. 


of the falling time of a stellar cluster. For a system of 
1.6 10® Mq, it would take a Hubble time to fall in the cen¬ 
tre from a radius of a few kpc. This timescale is strongly 
dependent on the initial radius r^, so that only ~ 10® yr 
are needed for the system to fall from a radius of 6.6 pc 
(with Vpi ~ 300 km.s”^, see Table ^ and Sect. 5.1). 

However, the application of the Chandrasekhar for¬ 
mula here is not appropriate. On one hand, this formula 
overestimates the dynamical friction since it assumes a 
singular isothermal sphere while PI has entered the core 
where the density remains constant. Moreover, the grav¬ 
itational potential is dominated by a central point-like 
mass, which does not provide friction. 


On the other hand, our simulations confirm a much 
more important result already noticed in the case of inter¬ 
acting galaxies by Prugniel & Combes (1992): the Chan¬ 
drasekhar formula underestimates the effect of dynami¬ 
cal friction since it assumes that the accreted system is 
rigid. Prugniel & Combes (1992) showed that releasing 
the hypothesis of a rigid satellite increases considerably 
the friction efficiency. In our experiments, the main dy¬ 
namical friction is indeed due to the deformation of the 
stellar cluster, a reciprocal effect usually neglected in es¬ 
timations based on the Chandrasekhar formula. 


run 

Mf,h 

g{MIL) 

rc 

Ro 

Vo 

#1 

6.5 10' 

1 

0.12 

00 

bo 

Vc 

#2 

9.0 10^ 

1 

0.12 

6.6 

Vc 

#3 

1.1 10® 

1 

0.12 

6.6 

Vc 

#4 

1.1 10® 

2 

0.12 

6.6 

Vc 

#5 

1.1 10® 

1 

0.12 

6.6 

1.5 X 14 

#6 

1.1 10® 

1 

0.12 

6.6 

Vc 

#7 

1.1 10® 

1 

0.12 

6.6 

-Vc 

#8 

1.1 10® 

1 

0.12 

6.6 

V. 

#9 

2.2 10^ 

1 

0.12 

6.6 

Vc 


4. The disruption of the cluster 

If PI indeed corresponds to an additional cold stellar sys¬ 
tem and is spatially close to the centre of the nucleus, its 
decaying time should be rather short, and the probabil¬ 
ity to observe such a configuration is very small. This is 
the main argument against the falling cluster hypothesis. 
However Lauer et al (1993) emphasize that the dynamical 
friction effects are not well known, and there exist some 
circumstances where it may shut off. They privileged the 
hypothesis of PI being an external stellar system coming 
from another nucleus, possibly having its own central BH. 

4.I. Dynamical friction 

A simple application of the Chandrasekhar formula (see 
e.g. Binney & Tremaine 1987) can give us a first estimate 


4-2. Tidal forces and the black hole mass 

Our simulations show that the main effect is due to the 
tidal forces induced by the central mass concentration 
which rapidly disrupt the cluster. The strength of the 
tidal forces strongly depends on the ratio of the central 
cluster and nucleus concentrations. For a black hole mass 
of ~ 6.10^ Mq, the cluster looses 30% of its mass in less 
than 10® yr. After 4.10® yr, they are all spread in a ring¬ 
like structure (Fig. ||), although the centre of gravity of the 
cluster’s particles is close to the centre of the nucleus. For 
masses greater than 6 10^ Mq, the black hole dominates 
the overall gravitational potential in the central 10 pc. The 
lifetime of the cluster is then roughly scaled like the inverse 
of the square root of the central dark mass (Fig. We 
also ran an experiment equivalent to run ffS but without 
the additional central mass. In this run, the stellar clus¬ 
ter is not disrupted. On the contrary, it is the nucleus 
now that suffers strong tidal deformations, and the clus¬ 
ter decays like a rigid body through dynamical friction, 
produced by the nucleus deformations. The decay time is 
then longer, of the order of 3 10® yr. 

In Fig. ^ we can also see that there is very little trans¬ 
fer of angular momentum between the nucleus and the 
original cluster: the mean angular momentum had only 
changed by 10% by the end of the simulation at 4.6 10® 
yr. This result is very similar to the one obtained by Charl¬ 
ton & Laguna (1995) at a larger scale: the debris of the 
globular clusters disrupted by tidal forces mainly follow 
the original cluster orbit (i.e. highly eccentric orbits). 
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Fig. 3. Face-on view of the cluster particles at different times 
of the simulation (run #1): the surface brightness has been 
slightly smoothed and only the faintest contours have been 
drawn to emphasize the development of the ring-like structure. 
The bottom right panel shows the histogram of the angular 
momentum for the cluster particles at the begining and at the 
end of the simulations (0.01 Myr - solid line - and 0.46 Myr - 
dashed line - respectively). 

In most runs presented in this paper, we have launched 
the cluster with the circular velocity. In run #5, we have 
launched the cluster with a tangential velocity of 1.5 x 14 
(hyperbolic orbit) and in run ^8 with a purely vertical ve¬ 
locity of 250 krn.s”^. In both cases, the distance between 
the cluster and the central mass concentration increases 
with time. This prevents the cluster to be strongly dis¬ 
rupted as illustrated in Fig. || (see also Charlton & Laguna 
1995). 

4-3. The mass-to-light ratio of the cluster 

Because of the absence of colour gradients, we assumed 
until now that PI and P2 have similar stellar populations 
and mass-to-light ratios. We have however also explored 
different mass-to-light ratios (hereafter M/L). For the run 
^4, we fixed the M/L of the cluster to be half the one of 



Fig. 4. Effect of the central mass concentration on the disrup¬ 
tion of the cluster: the mass of the cluster (in % of its total 
mass) enclosed in a sphere of radius r (normalized by the ini¬ 
tial core radius r^) is shown as the simulation evolves (from 
top to bottom, t = 10"^, 2.10"*, 3.10"^ and 4.10“*), for two differ¬ 
ent central mass concentrations (2.10^ 1^0“ solid lines - and 
10® Mq- dotted lines). 



Fig. 5. Same as Fig. ^ but for two different initial velocities of 
the cluster: run 4f5 corresponds to a cluster on an hyperbolic 
orbit in the equatorial plane, and run #8 to an initial pure 
vertical velocity for the cluster. Note the difference with run 
#3 in Fig. I 


the nucleus (see Table ^). This was achieved by assigning 
fewer particles to the cluster, which were then used for the 
nucleus whose total mass was conserved. The less massive 
cluster is obviously less concentrated at all times after 
the launch. However, its envelope seems to survive longer. 
To understand this point it is necessary to remember that 
the initial total energy of the clusters are different in those 
two runs: in run ^3 the cluster covers a longer track which 
induces a larger spreading of the unbound particles. More¬ 
over, the tidal radius of the cluster is significantly smaller 
(by about 25%) in run #4 than in run #3. 
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5. The case of M31 nucleus 


5.1. Some simple geometrieal arguments 


Let us first assume that the stellar cluster is located in 
the equatorial plane of the (axisymmetric) nucleus, and 
let (a:^, {V^, Vy) be its position and velocity in a refer¬ 

ence system centred on P2 where Ox is the line of nodes. 
Since we observe PI at about 0'.'14 of P2 along the nu¬ 
clear minor-axis and 0'.'5 along the major axis (Lauer et 
al. 1993, BEMN94), we can write: y‘^ = (y.'14/cosj with 
i being the inclination angle of the nucleus (90° = edge- 
on). In this case, the nucleus of M 31 cannot be viewed 
exactly edge-on. The average axis ratio of the nucleus 
of {q) of 0.61 can be used to bracket its inclination an¬ 
gle: arccos(0.61) ~ 52.4 < i{°) < 90. Then, the appar¬ 
ent distance of PI to the major-axis of the nucleus pro¬ 
vides us a lower limit for the true distance of PI to the 
line of nodes: ?/° > 0"14/0.61 ~ 0'.'23. And we have now 
-I- ?/°^) > 0"55 (~ 2 pc at 0.77 Mpc). If the 
near side of the nucleus corresponds to the near side of the 
main disc of M 31, i.e. the North-West side (see Fig. 3 of 
Kormendy 1988), we can conclude from the sign of the ob¬ 
served velocities that PI is going toward the line of nodes 
of the nucleus. 

The projection on the sky gives Vy = sini, where 
is the line-of-sight velocity of PI. If we assume now 
that the cluster is on a circular orbit, we can also write: 


W (R) 


X 


(3) 


with the spatial velocity of PI. Using Fig. 7 of Ko¬ 
rmendy & Richstone (data of Kormendy & Bender 1995), 
we estimate 1/)°^, to be ~ 160 km.s“^(see Sect. ||). This 
value is very uncertain, but the following calculations can 
be adapted accordingly. The projected distance x' (= x) 
along the nuclear major-axis is ~ 075 (Lauer et al. 1993, 
BEMN94). Therefore U° (i?) ~ 320-(y/((L25T^/sini). 
U° {R) is then minimum (U° ~ 211 km.s”^) for i ~ 62.6° 
which gives y^ ~ 073 and i?° ~ 07 59. 

In fact, since we already have observed values for a;° 
and we only need to get (or the inclination i) and 
14 ° to know the full position and velocity of the cluster. 
Using the model given in Bacon et al. (1994) for the spa¬ 
tial mass distribution and assuming a supermassive black 
hole of Mbh = 10® MQ(see Sect. ||), it is then possible to 
roughly constrain the type of orbit of the cluster in terms 
of these two unknowns and (as a test particle in 
a fixed potential). This is illustrated in Fig. where the 
limits of the regions for elliptic, and hyperbolic orbits are 
drawn as well as the curves for tangential and radial veloc¬ 
ities. We thus find that there is only one point where PI 
can be on a circular orbit: U° ~ 300 km.s~^for i?° ~ 1.05" 
and ~ 079 (see Fig. ^). Taking into account the uncer¬ 
tainty on the true value of on the circular velocity 

in the nucleus. Fig. ^ indicates that 0.85 < i?°(") < 1.3 or 
0.7 < y%”) < 1.2 (2.6 < y%pc) < 4.5 at 0.77 Mpc). 




Fig. 6. Diagrams showing the different types of orbits for a test 
particle in the modelled potential of M 31 (see text) in terms 
of i (or 7?°) and V^: the dashed lines correspond to the escape 
velocity, the dotted line to the circular velocity. Geometrical 
constraints are given by the solid line (pure tangential velocity) 
and the dotted-dashed line (radial orbits). Vertical dotted lines 
indicate the circular velocity solution. 


The core of the cluster could survive only if its cen¬ 
tral density is larger than the average mass density of the 
nucleus interior to its orbit. We have seen in the simula¬ 
tions that the initial central density of the cluster had to 
be « 10 times the average density at Rq = 6pc (includ¬ 
ing a black hole mass of 7 10^ M©). We therefore expect 
that the disruption will occur at a radius of ~ 3 pc. This 
roughly corresponds to the radius derived for a circular 
orbit inside the nucleus. If the stellar cluster follows such 
an orbit it is then presently being disrupted. Such quasi 
circular orbits were tried, where the cluster was slowly spi¬ 
raling in towards the centre because of dynamical friction. 
But the stripping of its outer parts was so effective that, 
although the center of mass of the system did reach the 
required radius, the remaining core of the globular clus¬ 
ter did not. Also, the disruption of the cluster produced a 
much too large internal velocity dispersion to the stellar 
debris around the core; the presently observed low veloc¬ 
ity dispersion of PI strongly suggests that this is not the 
case yet. Therefore, we must conclude that PI cannot be 
on a circular orbit around P2. 

Finally, for a central dark mass of 1. 10^ MQthe circu¬ 
lar velocity predicted from the model given in BFMN94 is 
smaller than 160 km.s~^at radii larger than 0755. There¬ 
fore, if the mass concentration is smaller than 1. 10^ AIq, 
the cluster is on an orbit unbound with respect to the nu¬ 
cleus, and will escape from its gravitational influence (and 
might eventually fall back later by friction). 

We have not explored many different inclinations of the 
orbital plane. However, as the main influence comes from 
the central dark mass, the fate of the cluster is mainly 
determined from its pericentre distance. 

5.2. Observations of the experiments 

The main problem to fit self-consistently the nuclear disc 
of M31 is its observed thickness, mainly under the assump- 












E. Emsellem & F. Combes: N body simulations of the nucleus of M 31 


tion of its nearly edge-on inclination. This requires a rel¬ 
atively high z-dispersion, while the disc is highly rotating 
in the plane. The system relaxes quickly, on time scales of 
the order of 10^ yr, and it revealed very difficult to main¬ 
tain both a rather thick disc, and a high V,.ot/cr ratio. 
With relaxation, the velocity dispersion tends to become 
more isotropic, and the rotational velocity decreases with 
respect to the in-plane dispersion. To reproduce the high 
azimutal velocities, a more massive black hole is required, 
and no good fit was obtained below 8 10 ^ Mq. 

We have “observed” the simulations by associating a 
specific luminosity of 50.5 LQper particle as to reproduce 
the photometric data in the / band (i.e. M/Lj ~ 2.5). 
For the bulge, we used the axisymmetric surface bright¬ 
ness model given in BEMN94 (including the difference in 
position angle, see their Table 4). The system (nucleus + 
cluster) was viewed such as the strongest projected peak 
of the cluster is located at x' = 1.85 pc and y' = 0.5 pc 
from the centre of the nucleus (see BEMN94). This fixes 
the three Euler angles of the projection. All the 

relevant characteristics of the various instrumental set¬ 
ups are then taken into account: spatial and spectral pixel 
sizes, spatial and spectral resolutions. The Point Spread 
Function was approximated by a single gaussian function 
which is thus completely fixed by its width tr*. For the 
bulge we assumed a gaussian Line Of Sight Velocity Dis¬ 
tribution (hereafter LOSVD) whose first two moments are 
those of the corresponding isotropic axisymmetric model 
(i.e. Model A of BEMN94 with M/Li = 2.5 and the ap¬ 
propriate central dark mass). The obtained total LOSVDs 
were parametrized with the help of Gauss-Hermite func¬ 
tions (see van der Marel & Franx 1993) 

Let us first examine the hypothesis of a projection ef¬ 
fect, where the cluster is physically far from the centre 
of M 31. The results of these experiments will be used to 
constrain a fully self-consistent simulation. 

5.3. Stellar cluster seen in projection against the nucleus 

Lauer et al. (1993) argued that the observed photometry 
and kinematics speaks against a cluster physically far from 
the centre and seen in projection: 

— Pi’s luminosity profile is truncated and asymmetric: 
this is indeed suggestive of an interaction between PI 
and P2. However, as mentioned by King et al. (1995), 
the decomposition is certainly not unique 

— Lauer et al. (1993) also argued that the positioning 
of PI on the nucleus major-axis favors a physical re¬ 
lationship between PI and P2. But this assumes that 
P2 is seen almost edge-on. Moreover, PI is closer to 
the major-axis of the bulge. 

— The velocity profiles are symmetric, as though PI did 
not introduce any perturbation. But this could come 
in part from the fact that the observed kinematics is 
luminosity weighted (seeing effect, cf Sect. |^ ). 


~ Lauer et al. (1993) estimated that the PI velocity dis¬ 
persion is too large for a globular cluster. However, the 
recent SIS data obtained by Kormendy indicate that 
api < 85 km.s“^(KR95). Moreover, our simulations 
suggest (see Sect. ^.3.1 ) that 30 < (Tpi(km.s“^) < 80 
(if PI is being disrupted it is certainly not in virial 
equilibrium). 

Although PI and P2 seem to be interacting, there are thus 
no definite arguments against a projection effect. We will 
therefore examine this possibility in the following para¬ 
graphs, and use this to test the stability of the model, 
after a relaxation of several dynamical times. 


The apparent flattening of the nucleus (see Sect. p.l[) 
sets a lower limit of ~ 53° for its inclination angle. We have 
attempted to reproduce the photometry and the kinemat¬ 
ics of M31 ’s nucleus using a thin disc and a central dark 
mass of a few 10^ Mq. We managed to get a reasonable 
fit with Mbh = 7-7 10^ MQand an inclination of 55°. The 
cluster has been fixed as a separate Plummer sphere with 
Tc = 0.8 pc, projected onto the nucleus. The simulation of 
the nucleus has been continued until the system reached 
dynamical stability. Self-consistent experiments are de¬ 
scribed in Sect. 

In Fig. 1^ we present the isophoteJl of the model com¬ 
pared to the WFPC2 / band image kindly provided 
by Michael Rich. We normalized both images using the 
HRCAM Ic band image (BEMN94). The fit is rather good, 
although we do not reproduce the flattening of PI as well 
as its full spatial extension (see King et al. 1995). Our nu¬ 
cleus is also slightly more peaked than in the observations. 
This is due to the fact that we start from conditions that 
are close to what is observed. Then dynamical relaxation 
is at work and the nucleus has a tendency to concentrate. 



Fig. 7. A band isophotes of the “inclined” model (left panel) 
and of the WFPC2 I band image (normalized to the Cousin 
system, right panel). The brightest isophote corresponds to 
11.85 mag.arcsec”^, and the step is 0.3 mag.arcsec”^. 


® In all the following figures, we have aligned the nucleus 
major-axis with the horizontal axis: see BEMN94. 
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5.3.1. Comparison with published kinematics 

As discussed above, the line-of-sight velocity of the clus¬ 
ter suggested by the available observables is probably 
close to 160 km.s”^. However, it is difficult to determine 
^los certainty. A way to constrain the mean value 

for is to directly look at the influence of PI on the 
LOSVDs. This has been already done by BEMN94 who 
unfortunately only examined a very limited subset of pos¬ 
sibilities. Moreover, at the spatial and spectral resolution 
of the TIGER data (cr* ~ 0"38) the light contribution of 
PI is significantly diluted. BEMN94 still observed a slight 
asymmetry in the velocity curve along the nucleus major- 
axis which hinted for a rather high apparent value of 
(> 120 km.s”^). We have fixed to a number of differ¬ 
ent values to test for its effect on the shape of the resulting 
LOSVDs and the corresponding derived velocities and dis¬ 
persions. The internal velocity dispersion of the cluster is 
assumed isotropic with cr'^ ~ 36 km.s”^. 

For comparison we have used the kinematics ob¬ 
tained with the two-dimensional TIGER data presented 
by BEMN94 (spatial resolution of cr* ~ 0'.'37 and spectral 
resolution of ~ 100 km.s”^). The best fit is obtained for 
VjQg = 140 km.s~^but values of 120 and 160 km.s~^are 
also consistent with the data (Fig. ||). Higher velocities 
for the cluster tend to increase both the peak in the ve¬ 
locity and velocity dispersion profiles, but to reduce the 
influence on higher order moments. It is important to 
remember that the kinematics were derived here from 
a gaussian fit to the LOSVDs, and that true statistical 
moments are less influenced by the presence of the clus¬ 
ter. Although Fig. H only shows the major-axis kinemat- 
ical profiles, the agreement is excellent on the full two- 
dimensional fields. The higher resolution data of Kor- 
mendy (SIS/GFHT, cr* = 0'.'27; KR95) also suggests a 
value of ~ 160 krn.s”^, which will therefore be used 
for the “self-consistent” model (Sect. ^a|) . 

We have finally estimated the effect of a higher ve¬ 
locity dispersion of the cluster by increasing its mass by 
a factor of 2: this increases the internal velocity disper¬ 
sion by a factor of v^. As expected the velocity pro¬ 
files are almost coincident but the velocity dispersion pro¬ 
files differ significantly: the minimum along the major-axis 
reaches ~ 95 km.s~^for tr'^ = 36 km.s~^(M/L/ = 2.5) and 
~ 112 km.s~^for cr“ = 51 km.s~^(M/L/ = 5) closer to the 
observed values. This suggests a reasonable interval for the 
internal dispersion of the cluster: 30 < cr'^(km.s“^) < 80. 

5.4- Stellar cluster inside the nucleus 

We have then examined the case of a rather more edge-on 
model for the nucleus, but this time including the stellar 
cluster in a self-consistent way. In this case, the nucleus 
should be significantly thicker as the apparent flattening 
is closer to the true one, and we have chosen a z veloc¬ 
ity dispersion slightly higher than the radial dispersion 





Fig. 8. TIGER velocity, dispersion, /is, and hi major-axis pro¬ 
files (circles). The corresponding profiles derived from the “in¬ 
clined model” are presented for different values of Vi^^: 100 and 
120 km.s“^(dashed lines), 140 km.s“^(solid line), 160 and 180 
km.s“^(dotted lines). 


{cTz = l.3<Tr, see Sect. 3.2.2). The nucleus is projected with 
an inclination ~ 70°. The orbit of the stellar cluster is no 
longer in the equatorial plane, but follows a nearly planar 
orbit, inclined by ^ 20°, so that the orbit is almost edge-on 
projected on the sky. We have already excluded a circu¬ 
lar orbit dO), and therefore have chosen an hyperbolic 
orbit, the cluster being now very close to its pericenter 
of 3pc, where it has barely reached the critical radius of 
disruption. 

We obtained a reasonable fit to the observables for a 
central dark mass of 9.4 10^ MQand an inclination of 68° 
(constrained by the observed position of the cluster). The 
cluster was launched at 8.8 pc from the centre with a ve¬ 
locity of 420 krn.s”^. The cluster has a strong effect on the 
stability of the nucleus: the outer parts gain energy and 
inflate, while the center becomes more condensed. This is 
difficult to anticipate in the initial conditions, and causes 
the main discrepancy between this model and the observ¬ 
ables: the central surface brightness of P2 is larger in the 
model than in the HST photometry (Fig. ^). Also, the 
nucleus has the tendency to dynamically relax towards a 
more centrally concentrated and flattened system, due to 
the strong influence of the central dark mass. The latter 
completely dominates the gravitational potential in the 
inner few arcseconds and tends to symmetrize the disper¬ 
sion tensor in the centre. We should note however that 
although P2 is too concentrated in the model it has the 
correct integrated magnitude. 
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Fig. 9. Same figure as Fig. but for the “self-consistent” 
model. 


Figs. 0011 and 1^ present the comparison of this 
model observed at t = 1.6 10^ years after the launch of 
the cluster with the TIGER and SIS data. The agreement 
is rather good even for the higher order Gauss-Hermite 
moments and hi (Fig. p^ ). As mentioned above, the 
isophotes are slightly too flattened which leads to a higher 
central surface brightness than in the observed image. At 





Fig. 10. Same as Fig. ^ but for the self-consistent” model. 


a time t = 1.6 10^ years, the cluster which lies at about 
1"1 (4.1 pc) is already partially disrupted, although the 
core still survives with a projected FWHM of ~ 0'.'4. This 
is slightly smaller than the value of 0'.'5 FWHM derived by 
King et al. (1995). The spread of the particles of the cluster 
induced by the tidal forces is accompanied by the creation 
of a large internal velocity gradient which in projection 
goes from ^ —105 km.s“^to ^ 240km.s“^. This dilutes 
the contribution of the cluster to the observed LOSVDs. 




Fig. 11. Isophotes (bottom), velocity fields (left panels) and 
velocity dispersion fields (right panels) from the TIGER data 
(top panels) and of the “self-consistent model” (middle panels). 
The steps are respectively 0.25 magnitude, 25 km.s“^(F) and 
15 km.s“^(cr). 


5.5. M 31 at the HST resolution 

At the HST resolution, the observed kinematics of the nu¬ 
cleus of M 31 should clearly reveal the presence of a super¬ 
imposed cold cluster. This is illustrated in Fig ^ where we 
have derived the kinematics of the “self-consistent model” 
assuming a pixel size of 0'.'06 x 0'.'06. We also present 
seven corresponding LOSVDs_along the line joining PI 
and P2 for squared aperturesQ of 0'.'26 and (y.'09 and a 
spectral resolution of 100 km.s”^. Note the large wings of 
the LOSVD at the central position due to the presence 
of the mass concentration. There are no significant dif¬ 
ferences between the observed LOSVDs for the two differ¬ 
ent apertures. The velocity (gaussian fit) along the major- 

^ These parameters were chosen to be similar to the charac¬ 
teristics of the observations of M 31 with the HST Faint Object 
Spectrograph. 
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Fig. 12. Velocity and velocity dispersion of the SIS data (cir¬ 
cles, KR95 and of the “self-consistent model”. 



Fig. 13. LOSVDs derived from the “self-consistent model” at 
four different positions to be compared with Fig. 7 of KR95. 

axis reaches ~ —220 km.s~^on the side opposite to PI at 
~ 0'.'45, and the maximum dispersion is ~ 290 km.s~^at 
~ 0V25. The maximum positive velocity and the minimum 
dispersion obviously strongly depend on the true internal 
kinematics of PI. The main prediction for spectrographic 
observations at the HST resolution is the strong asymme¬ 
try in the velocity profile along the P1-P2 line, which is 
also clearly revealed in the shape of the LOSVDs. How¬ 
ever, whatever PI corresponds to (a stellar cluster or a 
density perturbation), the fact that the observed kine¬ 
matics is weighted by light induce strong asymmetries. 
Finally, as noticed by Lauer et al. (1993), if the cluster 
has indeed a transverse velocity of ~ 440 km.s~^as in the 
“self-consistent model” presented here, its proper motion 
should be ~ 10“^ arcsecond in 10 years. 

6. Discussion and conclusions 

In this paper we have presented N-body simulations of a 
stellar cluster falling towards a galactic nucleus. These ex¬ 
periments were designed to ressemble the central region of 
M 31. First we have shown that it was possible to obtain 
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Fig. 14. LOSVDs from the “self-consistent model” at —Cf.'75, 
-0'.'5, -0'.'25, 0, -f0'.'25, -0'.'5, -f0'.'75 along the P1/P2 line for 
squared apertures of 0”26 (solid lines) and 0”09 (dotted lines). 




Fig. 15. Velocity (left) and dispersion profiles (right) along the 
major-axis of the nucleus (solid lines) and along the P1/P2 line 
(dotted lines) predicted from the “self-consistent model”. 

reasonable fits to all observables, either with a projected 
stellar cluster (which served to adapt our initial conditions 
for the N body simulations), or with a stellar cluster close 
to the nucleus and interacting with it. The best model 
considered an initial orbit where the cluster was unbound 
to the nucleus. Slight discrepancies between this model 
and the data still remain (e.g. the nucleus isophotes flat¬ 
tening). This solution is certainly not unique and other 
significantly different initial conditions (e.g. the choice of 
the orbit) could lead to very similar observables. 

Second, we have shown that the dynamical friction is 
even more efficient than previously estimated: the cluster 
rapidly decays to smaller distances, where it is disrupted 
through tidal interaction with the BH. The lifetime of such 
a cluster is short (< 5 10^) at radii smaller than ~ 10 pc if 
a central dark mass of a few 10^ MqIs present at the cen¬ 
tre. The main argument against the scenario of a falling 
cluster is that we must see the M31 nucleus at a very spe¬ 
cial time. However, the decaying of a stellar cluster is not 
a rare phenomenon. Tremaine et al (1975) proposed that 
the falling of globular clusters could even be the formation 
mechanism of the nuclear disc. Orders of magnitude are 
consistent with the observations, if about 25 clusters have 
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been disrupted to form the nucleus, which mass has grown 
as This means that one cluster is falling every 400 
Myr; since it will perturb the brightness distribution dur¬ 
ing « 0.2 Myr, the probability of the present configuration 
of the M31 nucleus is about 4 10“^. This however includes 
only globular clusters spiraling in on a quasi-circular or¬ 
bit. We have shown that PI being on a circular orbit is 
inconsistent with the presence of a central dark mass of a 
few 10^ Mq. More rapid globular clusters, in hyperbolic 
orbits, as in one of the scenarii proposed here, would be 
disrupted, but not slowed down enough to contribute to 
the nucleus mass. 

Therefore, we should conclude that the hypothesis of 
a falling cluster (e.g. nucleus of an accreted dwarf galaxy, 
globular cluster) remains a viable explanation. Lauer et 
al. (1993) suggested that PI could survive close to the 
centre if it contains a secondary black hole of “substantial 
mass”. The probability of seeing it now in this configu¬ 
ration is even smaller than for a stellar cluster, given the 
frequencies of black holes. Only if we assume that dark 
matter in galactic halos is constituted by BHs are the lat¬ 
ter more numerous than globular clusters. In the model 
of Lacey & Ostriker (1985), one of these 10® Mq BH ar¬ 
rives in the centre every 10® yrs, and most of the time it 
merges with the BH already present, or is ejected if there 
is already a binary BH. Xu & Ostriker (1994) argue that, 
since fresh BH arrive regularly on the centre, the proba¬ 
bility to find a binary BH with a separation of 1.6pc (as 
P1-P2) is 10%. It is however difficult to explain why PI is 
so bright, without any cusp contrary to P2. More stringent 
is the constraint of a low dispersion at PI: in the case of 
a secondary black hole at PI, the dispersion there should 
be signihcantly larger than observed. 

Our simulations hnally served to precise the required 
BH mass, in the hypothesis of an axisymmetric nucleus: in 
order to ht the disc thickness and stability requirements, 
we need the presence of a mass concentration of at least 
7 10^ MQto explain the central stellar velocity and veloc¬ 
ity dispersion gradients. The detection of a strong central 
UV and X ray source in the centre of M 31 is also in favor 
of a supermassive black hole. The main underlying as¬ 
sumption in the kinematical determination of the central 
dark mass is the axisymmetry of the nucleus. In a forth¬ 
coming paper, we will examine the possibility of a triaxial 
morphology (i.e. a bar) for the nucleus of M 31. 
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